// nvcc -g -lglut -lGLEW  -I$CUDA_SDK/C/common/inc -L$CUDA_SDK/C/lib -lcutil_i386 whitenoise.cu
// Compile ^^^
// ========================================================================
// // Kurs: Procesory graficzne w obliczeniach równoległych 
//          A.Łukaszewski 2010
//=========================================================================
// CUDA-OpenGL interoperability ===========================================
// based on  CUDA SDK : SimpleGL             PBO - pixel buffer object here  
//=========================================================================
#include <stdio.h>
#include <GL/glew.h>
#include <GL/glut.h>
#include <cutil.h>
#include <cuda_gl_interop.h>
#include <string>
#include <iostream>

__device__ uint d_width = 800, d_height = 600;
           uint   width = 800,   height = 600;

typedef unsigned int  uint;

void reset_vars();

// http://fly.cc.fer.hr/~unreal/theredbook/chapter08.html
GLubyte rasters[][13] = {
{0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0x18, 0x18, 0x00, 0x00, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18}, 
{0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x36, 0x36, 0x36, 0x36}, 
{0x00, 0x00, 0x00, 0x66, 0x66, 0xff, 0x66, 0x66, 0xff, 0x66, 0x66, 0x00, 0x00}, 
{0x00, 0x00, 0x18, 0x7e, 0xff, 0x1b, 0x1f, 0x7e, 0xf8, 0xd8, 0xff, 0x7e, 0x18}, 
{0x00, 0x00, 0x0e, 0x1b, 0xdb, 0x6e, 0x30, 0x18, 0x0c, 0x76, 0xdb, 0xd8, 0x70}, 
{0x00, 0x00, 0x7f, 0xc6, 0xcf, 0xd8, 0x70, 0x70, 0xd8, 0xcc, 0xcc, 0x6c, 0x38}, 
{0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x18, 0x1c, 0x0c, 0x0e}, 
{0x00, 0x00, 0x0c, 0x18, 0x30, 0x30, 0x30, 0x30, 0x30, 0x30, 0x30, 0x18, 0x0c}, 
{0x00, 0x00, 0x30, 0x18, 0x0c, 0x0c, 0x0c, 0x0c, 0x0c, 0x0c, 0x0c, 0x18, 0x30}, 
{0x00, 0x00, 0x00, 0x00, 0x99, 0x5a, 0x3c, 0xff, 0x3c, 0x5a, 0x99, 0x00, 0x00}, 
{0x00, 0x00, 0x00, 0x18, 0x18, 0x18, 0xff, 0xff, 0x18, 0x18, 0x18, 0x00, 0x00}, 
{0x00, 0x00, 0x30, 0x18, 0x1c, 0x1c, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xff, 0xff, 0x00, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0x00, 0x38, 0x38, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x60, 0x60, 0x30, 0x30, 0x18, 0x18, 0x0c, 0x0c, 0x06, 0x06, 0x03, 0x03}, 
{0x00, 0x00, 0x3c, 0x66, 0xc3, 0xe3, 0xf3, 0xdb, 0xcf, 0xc7, 0xc3, 0x66, 0x3c}, 
{0x00, 0x00, 0x7e, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x78, 0x38, 0x18}, 
{0x00, 0x00, 0xff, 0xc0, 0xc0, 0x60, 0x30, 0x18, 0x0c, 0x06, 0x03, 0xe7, 0x7e}, 
{0x00, 0x00, 0x7e, 0xe7, 0x03, 0x03, 0x07, 0x7e, 0x07, 0x03, 0x03, 0xe7, 0x7e}, 
{0x00, 0x00, 0x0c, 0x0c, 0x0c, 0x0c, 0x0c, 0xff, 0xcc, 0x6c, 0x3c, 0x1c, 0x0c}, 
{0x00, 0x00, 0x7e, 0xe7, 0x03, 0x03, 0x07, 0xfe, 0xc0, 0xc0, 0xc0, 0xc0, 0xff}, 
{0x00, 0x00, 0x7e, 0xe7, 0xc3, 0xc3, 0xc7, 0xfe, 0xc0, 0xc0, 0xc0, 0xe7, 0x7e}, 
{0x00, 0x00, 0x30, 0x30, 0x30, 0x30, 0x18, 0x0c, 0x06, 0x03, 0x03, 0x03, 0xff}, 
{0x00, 0x00, 0x7e, 0xe7, 0xc3, 0xc3, 0xe7, 0x7e, 0xe7, 0xc3, 0xc3, 0xe7, 0x7e}, 
{0x00, 0x00, 0x7e, 0xe7, 0x03, 0x03, 0x03, 0x7f, 0xe7, 0xc3, 0xc3, 0xe7, 0x7e}, 
{0x00, 0x00, 0x00, 0x38, 0x38, 0x00, 0x00, 0x38, 0x38, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0x30, 0x18, 0x1c, 0x1c, 0x00, 0x00, 0x1c, 0x1c, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0x06, 0x0c, 0x18, 0x30, 0x60, 0xc0, 0x60, 0x30, 0x18, 0x0c, 0x06}, 
{0x00, 0x00, 0x00, 0x00, 0xff, 0xff, 0x00, 0xff, 0xff, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0x60, 0x30, 0x18, 0x0c, 0x06, 0x03, 0x06, 0x0c, 0x18, 0x30, 0x60}, 
{0x00, 0x00, 0x18, 0x00, 0x00, 0x18, 0x18, 0x0c, 0x06, 0x03, 0xc3, 0xc3, 0x7e}, 
{0x00, 0x00, 0x3f, 0x60, 0xcf, 0xdb, 0xd3, 0xdd, 0xc3, 0x7e, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0xc3, 0xc3, 0xc3, 0xc3, 0xff, 0xc3, 0xc3, 0xc3, 0x66, 0x3c, 0x18}, 
{0x00, 0x00, 0xfe, 0xc7, 0xc3, 0xc3, 0xc7, 0xfe, 0xc7, 0xc3, 0xc3, 0xc7, 0xfe}, 
{0x00, 0x00, 0x7e, 0xe7, 0xc0, 0xc0, 0xc0, 0xc0, 0xc0, 0xc0, 0xc0, 0xe7, 0x7e}, 
{0x00, 0x00, 0xfc, 0xce, 0xc7, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0xc7, 0xce, 0xfc}, 
{0x00, 0x00, 0xff, 0xc0, 0xc0, 0xc0, 0xc0, 0xfc, 0xc0, 0xc0, 0xc0, 0xc0, 0xff}, 
{0x00, 0x00, 0xc0, 0xc0, 0xc0, 0xc0, 0xc0, 0xc0, 0xfc, 0xc0, 0xc0, 0xc0, 0xff}, 
{0x00, 0x00, 0x7e, 0xe7, 0xc3, 0xc3, 0xcf, 0xc0, 0xc0, 0xc0, 0xc0, 0xe7, 0x7e}, 
{0x00, 0x00, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0xff, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3}, 
{0x00, 0x00, 0x7e, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x7e}, 
{0x00, 0x00, 0x7c, 0xee, 0xc6, 0x06, 0x06, 0x06, 0x06, 0x06, 0x06, 0x06, 0x06}, 
{0x00, 0x00, 0xc3, 0xc6, 0xcc, 0xd8, 0xf0, 0xe0, 0xf0, 0xd8, 0xcc, 0xc6, 0xc3}, 
{0x00, 0x00, 0xff, 0xc0, 0xc0, 0xc0, 0xc0, 0xc0, 0xc0, 0xc0, 0xc0, 0xc0, 0xc0}, 
{0x00, 0x00, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0xdb, 0xff, 0xff, 0xe7, 0xc3}, 
{0x00, 0x00, 0xc7, 0xc7, 0xcf, 0xcf, 0xdf, 0xdb, 0xfb, 0xf3, 0xf3, 0xe3, 0xe3}, 
{0x00, 0x00, 0x7e, 0xe7, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0xe7, 0x7e}, 
{0x00, 0x00, 0xc0, 0xc0, 0xc0, 0xc0, 0xc0, 0xfe, 0xc7, 0xc3, 0xc3, 0xc7, 0xfe}, 
{0x00, 0x00, 0x3f, 0x6e, 0xdf, 0xdb, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0x66, 0x3c}, 
{0x00, 0x00, 0xc3, 0xc6, 0xcc, 0xd8, 0xf0, 0xfe, 0xc7, 0xc3, 0xc3, 0xc7, 0xfe}, 
{0x00, 0x00, 0x7e, 0xe7, 0x03, 0x03, 0x07, 0x7e, 0xe0, 0xc0, 0xc0, 0xe7, 0x7e}, 
{0x00, 0x00, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0xff}, 
{0x00, 0x00, 0x7e, 0xe7, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3}, 
{0x00, 0x00, 0x18, 0x3c, 0x3c, 0x66, 0x66, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3}, 
{0x00, 0x00, 0xc3, 0xe7, 0xff, 0xff, 0xdb, 0xdb, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3}, 
{0x00, 0x00, 0xc3, 0x66, 0x66, 0x3c, 0x3c, 0x18, 0x3c, 0x3c, 0x66, 0x66, 0xc3}, 
{0x00, 0x00, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x3c, 0x3c, 0x66, 0x66, 0xc3}, 
{0x00, 0x00, 0xff, 0xc0, 0xc0, 0x60, 0x30, 0x7e, 0x0c, 0x06, 0x03, 0x03, 0xff}, 
{0x00, 0x00, 0x3c, 0x30, 0x30, 0x30, 0x30, 0x30, 0x30, 0x30, 0x30, 0x30, 0x3c}, 
{0x00, 0x03, 0x03, 0x06, 0x06, 0x0c, 0x0c, 0x18, 0x18, 0x30, 0x30, 0x60, 0x60}, 
{0x00, 0x00, 0x3c, 0x0c, 0x0c, 0x0c, 0x0c, 0x0c, 0x0c, 0x0c, 0x0c, 0x0c, 0x3c}, 
{0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xc3, 0x66, 0x3c, 0x18}, 
{0xff, 0xff, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x18, 0x38, 0x30, 0x70}, 
{0x00, 0x00, 0x7f, 0xc3, 0xc3, 0x7f, 0x03, 0xc3, 0x7e, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0xfe, 0xc3, 0xc3, 0xc3, 0xc3, 0xfe, 0xc0, 0xc0, 0xc0, 0xc0, 0xc0}, 
{0x00, 0x00, 0x7e, 0xc3, 0xc0, 0xc0, 0xc0, 0xc3, 0x7e, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0x7f, 0xc3, 0xc3, 0xc3, 0xc3, 0x7f, 0x03, 0x03, 0x03, 0x03, 0x03}, 
{0x00, 0x00, 0x7f, 0xc0, 0xc0, 0xfe, 0xc3, 0xc3, 0x7e, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0x30, 0x30, 0x30, 0x30, 0x30, 0xfc, 0x30, 0x30, 0x30, 0x33, 0x1e}, 
{0x7e, 0xc3, 0x03, 0x03, 0x7f, 0xc3, 0xc3, 0xc3, 0x7e, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0xc3, 0xfe, 0xc0, 0xc0, 0xc0, 0xc0}, 
{0x00, 0x00, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x00, 0x00, 0x18, 0x00}, 
{0x38, 0x6c, 0x0c, 0x0c, 0x0c, 0x0c, 0x0c, 0x0c, 0x0c, 0x00, 0x00, 0x0c, 0x00}, 
{0x00, 0x00, 0xc6, 0xcc, 0xf8, 0xf0, 0xd8, 0xcc, 0xc6, 0xc0, 0xc0, 0xc0, 0xc0}, 
{0x00, 0x00, 0x7e, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x78}, 
{0x00, 0x00, 0xdb, 0xdb, 0xdb, 0xdb, 0xdb, 0xdb, 0xfe, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0xc6, 0xc6, 0xc6, 0xc6, 0xc6, 0xc6, 0xfc, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0x7c, 0xc6, 0xc6, 0xc6, 0xc6, 0xc6, 0x7c, 0x00, 0x00, 0x00, 0x00}, 
{0xc0, 0xc0, 0xc0, 0xfe, 0xc3, 0xc3, 0xc3, 0xc3, 0xfe, 0x00, 0x00, 0x00, 0x00}, 
{0x03, 0x03, 0x03, 0x7f, 0xc3, 0xc3, 0xc3, 0xc3, 0x7f, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0xc0, 0xc0, 0xc0, 0xc0, 0xc0, 0xe0, 0xfe, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0xfe, 0x03, 0x03, 0x7e, 0xc0, 0xc0, 0x7f, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0x1c, 0x36, 0x30, 0x30, 0x30, 0x30, 0xfc, 0x30, 0x30, 0x30, 0x00}, 
{0x00, 0x00, 0x7e, 0xc6, 0xc6, 0xc6, 0xc6, 0xc6, 0xc6, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0x18, 0x3c, 0x3c, 0x66, 0x66, 0xc3, 0xc3, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0xc3, 0xe7, 0xff, 0xdb, 0xc3, 0xc3, 0xc3, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0xc3, 0x66, 0x3c, 0x18, 0x3c, 0x66, 0xc3, 0x00, 0x00, 0x00, 0x00}, 
{0xc0, 0x60, 0x60, 0x30, 0x18, 0x3c, 0x66, 0x66, 0xc3, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0xff, 0x60, 0x30, 0x18, 0x0c, 0x06, 0xff, 0x00, 0x00, 0x00, 0x00}, 
{0x00, 0x00, 0x0f, 0x18, 0x18, 0x18, 0x38, 0xf0, 0x38, 0x18, 0x18, 0x18, 0x0f}, 
{0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18}, 
{0x00, 0x00, 0xf0, 0x18, 0x18, 0x18, 0x1c, 0x0f, 0x1c, 0x18, 0x18, 0x18, 0xf0}, 
{0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x06, 0x8f, 0xf1, 0x60, 0x00, 0x00, 0x00} 
};

GLuint fontOffset;
void makeRasterFont(void) {
    GLuint i;
    glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
    fontOffset = glGenLists (128);
    for (i = 32; i < 127; i++) {
        glNewList(i+fontOffset, GL_COMPILE);
            glBitmap(8, 13, 0.0, 2.0, 10.0, 0.0, rasters[i-32]);
        glEndList();
    }
}

void printString(char *s) {
    glPushAttrib (GL_LIST_BIT);
    glListBase(fontOffset);
    glCallLists(strlen(s), GL_UNSIGNED_BYTE, (GLubyte *) s);
    glPopAttrib ();
}


const char a0 = 0;
const char a1 = 1;
const char a2 = 2;
const char a3 = 3;
const char b0 = 4;
const char b1 = 5;
const char b2 = 6;
const char b3 = 7;
const char n0 = 8;
const char fw = 9;
const char fh = 10;

float h_vars[20];
__device__ float *d_vars;
float keyb_step = 0.1;


//=========================================================================
// Pseudo Random: TausStep, LCGStep, Hybrid based on:
// GPU Gems 3:
// Lee Howes, David Thomas (Imperial College London)
// Chapter 37. Efficient Random Number Generation and Application Using CUDA
//=========================================================================
// Cheap pseudo random numbers:
//  
// S1, S2, S3, M - constants,  z - state
__device__ uint TausStep(uint &z, int S1, int S2, int S3, uint M)  {
    uint b=(((z << S1) ^ z) >> S2);
    return z = (((z & M) << S3) ^ b);
}

// A, C - constants
__device__ uint LCGStep(uint &z, uint A, uint C) {
    return z=(A*z+C);
}

// Mixed :
__device__ float HybridTaus(uint &z1, uint &z2, uint &z3, uint &z4) {
    // Combined period is lcm(p1,p2,p3,p4)~ 2^121
    return 2.3283064365387e-10 * (              // Periods
               TausStep(z1, 13, 19, 12, 4294967294UL) ^   // p1=2^31-1
               TausStep(z2,  2, 25,  4, 4294967288UL) ^   // p2=2^30-1
               TausStep(z3,  3, 11, 17, 4294967280UL) ^   // p3=2^28-1
               LCGStep( z4,    1664525, 1013904223UL)     // p4=2^32
           );
}

// Int Mixed and modified: cheaper
__device__ uint HybridTausInt(uint &z1, uint &z2, uint &z3, uint &z4) {
    // Combined period is lcm(p1,p2,p3,p4)~ 2^121
    return (              // Periods
               TausStep(z1, 13, 19, 12, 4294967294UL) ^   // p1=2^31-1
               //  TausStep(z2,  2, 25,  4, 4294967288UL) ^   // p2=2^30-1
               //  TausStep(z3,  3, 11, 17, 4294967280UL) ^   // p3=2^28-1
               LCGStep( z4,    1664525, 1013904223UL)     // p4=2^32
           );
}

// Testing func:   cheap one int state
__device__ uint funct(uint id) {
    //return LCGStep( id,    1664525, 1013904223UL) ;    // p4=2^32
    return HybridTausInt(id,id,id,id);
    //return id = (1664525*id+1013904223UL) % (65536*256);
    //return id = (xx%256) + 256*(y%256) + 65536*( (256-(xx%256)-(y%256))%256 ) ;
}
//=========================================================================

//=========================================================================
//initialization kernel:
__global__ void initim1(uint * d_output, uint imageW, uint imageH) {
    uint x  = __umul24(blockIdx.x, blockDim.x) + threadIdx.x;
    uint y  = __umul24(blockIdx.y, blockDim.y) + threadIdx.y;
    uint id = __umul24(y, imageW) + x;  // Good for < 16MPix

    if ( x < imageW && y < imageH ) {
        d_output[id] = id;
    }
}


//next pseudo random number:
__global__ void modify1(uint * d_output, uint imageW, uint imageH) {
    uint x  = __umul24(blockIdx.x, blockDim.x) + threadIdx.x;
    uint y  = __umul24(blockIdx.y, blockDim.y) + threadIdx.y;
    uint id = __umul24(y, imageW) + x;  // Good for < 16MPix

    // if ( x < imageW && y < imageH ) {
    // 	afin1_ntimes(&x, &y, 2);
    // 	d_output[id] = x*y;
    // }

    if ( x < imageW && y < imageH ) {
        d_output[id] = funct(d_output[id]);
	d_output[id] = 0; // start color
    }
}


/*
 * dla pary (x,y) liczy n razy (x,y)=f(a,b);
 */
__device__ void afin1_ntimes(uint *x, uint *y, uint n) {
    float newx;
    float newy;


    newx = float(*x) / d_width;
    newy = float(*y) / d_height;

    while(n--) {
	// uint random_value = newx*(*x);

	// if (random_value % 2) {
	//     // newx = -0.4f * (*x) - 1.0f;
	//     // newy = -0.4f * (*y) + 0.1f;
	//     newx = d_vars[a0] * (*x) + d_vars[a1];
	//     newy = d_vars[a2] * (*y) + d_vars[a3];
	// } else {

	//     // newx = 0.76f * (*x) -  0.4f * (*y);
	//     // newy = 0.4f  * (*x) + 0.76f * (*y);
	//     newx = d_vars[b0] * (*x) + d_vars[b1] * (*y);
	//     newy = d_vars[b2] * (*x) + d_vars[b3] * (*y);
	// }

	// newx = d_vars[a0] * (*x) + d_vars[a1];
	// newy = d_vars[a2] * (*y) + d_vars[a3];

	// newx = d_vars[a0] * newx + d_vars[a1];
	// newy = d_vars[a2] * newy + d_vars[a3];

	// newx = 0.76f * newx -  0.4f * newy;
	// newy = 0.4f  * newx + 0.76f * newy;

	    newx = d_vars[b0] * newx + d_vars[b1] * newy;
	    newy = d_vars[b2] * newx + d_vars[b3] * newy;


// 	newx = (*x) + 5;
// 	newy = (*y);
// 	}

    }
    *x = (uint)newx * d_width;
    *y = (uint)newy * d_height;
}


//next pseudo random number:
__global__ void ifs(uint * d_output, uint imageW,   uint imageH, 
		    float fractalW, float fractalH, uint iters, 
		    float ma0,      float ma1,     float ma2,   float ma3,
		    float mb0,      float mb1,     float mb2,   float mb3
) {
    uint x  = __umul24(blockIdx.x, blockDim.x) + threadIdx.x;
    uint y  = __umul24(blockIdx.y, blockDim.y) + threadIdx.y;
    uint id = __umul24(y, imageW) + x;  // Good for < 16MPix

//    uint color = d_output[id];
//
    
    float newx = float(x) / (float)d_width;
    float newy = float(y) / (float)d_height;
    float oldx = newx;
    float oldy = newy;
    
    uint random_value = id;
    uint start_point = id;
    uint n = iters;
    for(int i=0;i<n;i++){
	// random_value = funct(random_value);
	random_value = TausStep(start_point, 3,3,17,0xcafe);

	if (random_value % 8 < 5) {
	    // newx = -0.4f * oldx - 1.0f;
	    // newy = -0.4f * oldy + 0.1f;

	    // newx = d_vars[a0] * (*x) + d_vars[a1];
	    // newy = d_vars[a2] * (*y) + d_vars[a3];

	    newx = ma0 * oldx + ma1;
	    newy = ma2 * oldy + ma3;
	} else {

	    // newx = 0.76f * oldx -  0.4f * oldy;
	    // newy = 0.4f  * oldx + 0.76f * oldy;
	    newx = mb0 * oldx + mb1 * oldy;
	    newy = mb2 * oldx + mb3 * oldy;
	}

	oldx = newx;
	oldy = newy;
	x = (uint)((fractalW*newx+1.0f)*.5f * float(d_width));
	y = (uint)((fractalH*newy+1.0f)*.5f * float(d_height));


	// start
	if ( i>30 && x < imageW && y < imageH ) {
	    uint newid = __umul24(y, imageW) + x;  // Good for < 16MPix
	    d_output[newid] += 50;  // new color
	}

    }




    // if ( x < imageW && y < imageH ) {
    //     d_output[id] = funct(d_output[id]);
    // }

}


//=========================================================================
// Pseudo Random Kernels END
//=========================================================================


GLuint   pbo = 0;      // OpenGL PBO id.
uint    *d_output;     // CUDA device pointer to PBO data

dim3 blockSize(16,16); // threads
dim3 gridSize;         // set up in initPixelBuffer

int iDivUp(int a, int b) {
    return (a % b != 0) ? (a / b + 1) : (a / b);
}

void initPixelBuffer() {
    if (pbo) {      // delete old buffer
        cudaGLUnregisterBufferObject(pbo);
        glDeleteBuffersARB(1, &pbo);
    }
    // create pixel buffer object for display
    glGenBuffersARB(1, &pbo);
    glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, pbo);
    glBufferDataARB(GL_PIXEL_UNPACK_BUFFER_ARB, width*height*sizeof(GLubyte)*4, 0, GL_STREAM_DRAW_ARB);
    cudaGLRegisterBufferObject(pbo);

    // calculate new grid size
    gridSize = dim3(iDivUp(width, blockSize.x), iDivUp(height, blockSize.y));

    // from display:
    cudaGLMapBufferObject((void**)&d_output, pbo);
    initim1<<<gridSize, blockSize>>>(d_output, width, height);

    // losowe wartości
    modify1<<<gridSize, blockSize>>>(d_output, width, height);
    modify1<<<gridSize, blockSize>>>(d_output, width, height);

    CUT_CHECK_ERROR("Kernel error");
    cudaGLUnmapBufferObject(pbo);
}

static int cnt=0; // generation(display calls) count


__global__ void run_TausStep(uint * d_output, uint imageW,   uint imageH) {
    uint x  = __umul24(blockIdx.x, blockDim.x) + threadIdx.x;
    uint y  = __umul24(blockIdx.y, blockDim.y) + threadIdx.y;
    uint id = __umul24(y, imageW) + x;  // Good for < 16MPix

    if ( x < imageW && y < imageH ) {
        d_output[id] = funct(d_output[id]);
    }
} 

void check_generators(uint * d_output, uint imageW,   uint imageH) {
    // std::cout << "Checking generators" << std::endl;
    // modify1<<<gridSize, blockSize>>>(d_output, width, height);
    run_TausStep<<<gridSize, blockSize>>>(d_output, width, height);
    int size = imageW*imageH;
    int memsize = size*4;
    uint* h_taus = (uint*)malloc(memsize);
    cudaMemcpy(h_taus, d_output, memsize, cudaMemcpyDeviceToHost);
    free(h_taus);
}

void display() {

    check_generators(d_output, width, height);

    initPixelBuffer();

    // printf("%4d\n", cnt % 10000);
    cnt++;
    
    const int vars_size = 20 * sizeof(float);
    cudaMalloc(&d_vars, vars_size);
    cudaMemcpy(d_vars, h_vars, vars_size, cudaMemcpyHostToDevice);

    cudaGLMapBufferObject((void**)&d_output, pbo);
    ifs<<<gridSize, blockSize>>>(d_output, width, height, 
				 h_vars[fw], h_vars[fh], h_vars[n0],
				 h_vars[a0], h_vars[a1], h_vars[a2], h_vars[a3],
				 h_vars[b0], h_vars[b1], h_vars[b2], h_vars[b3]
	);
    CUT_CHECK_ERROR("Kernel error");
    cudaGLUnmapBufferObject(pbo );
    cudaThreadSynchronize();
    cudaFree(d_vars);

    glClear(GL_COLOR_BUFFER_BIT);
    glDisable(GL_DEPTH_TEST);
    glRasterPos2f(0.0f, 0.0f);
    glDrawPixels(width, height, GL_RGBA, GL_UNSIGNED_BYTE, 0);

    GLfloat white[3] = { 1,1,1 };
    glColor3fv(white);

    char astr[30]; 
    // przekształcenie (a) po lewej stronie
    {
	glRasterPos2f(0.01, 0.05);
	sprintf(astr, "x' = %+4.3fx %+3.3f", h_vars[a0], h_vars[a1]);
	printString(astr);
    }
    {
	glRasterPos2f(0.01, 0.01);
	sprintf(astr, "y' = %+4.3fy %+3.3f", h_vars[a2], h_vars[a3]);
	printString(astr);
    }

    // przeksztalcenie (b) po prawej stronie
    {
	glRasterPos2f(0.41, 0.05);
	sprintf(astr, "x' = %+4.3fx %+3.3fy", h_vars[b0], h_vars[b1]);
	printString(astr);
    }
    {
	glRasterPos2f(0.41, 0.01);
	sprintf(astr, "y' = %+4.3fx %+3.3fy", h_vars[b2], h_vars[b3]);
	printString(astr);
    }

    // zmienna n0 - liczba pominiętych elementów
    {
	glRasterPos2f(0.01, 0.01 + 0.05 * 3);
	sprintf(astr, "  iters   = %+3.0f", h_vars[n0]);
	printString(astr);
    }

    // zmienna keyb_step - krok zmiany zmiennych
    {
	glRasterPos2f(0.01, 0.01 + 0.05 * 2);
	sprintf(astr, "keyb_step = %+3.3f", keyb_step);
	printString(astr);
    }

    glFlush ();
    glutSwapBuffers();
    glutReportErrors();
}

void keyboard(unsigned char k, int , int ) {
    initPixelBuffer();
    if (k==27 || k=='Q') exit(1);
    glutPostRedisplay();

#define MANIP(var, keyincr, keydecr) if (k==keyincr) h_vars[var]+=keyb_step; \
                                     if (k==keydecr) h_vars[var]-=keyb_step
    MANIP(a0, 'q', 'w');
    MANIP(a1, 'e', 'r');
    MANIP(a2, 'a', 's');
    MANIP(a3, 'd', 'f');

    MANIP(b0, 't', 'y');
    MANIP(b1, 'u', 'i');
    MANIP(b2, 'g', 'h');
    MANIP(b3, 'j', 'k');
    MANIP(fw, '1', '2');
    MANIP(fh, '1', '2');
#undef MANIP

    if (k=='z') keyb_step+=0.005f;
    if (k=='x') keyb_step-=0.005f;

    if (k=='n') h_vars[n0]+=3;
    if (k=='m') h_vars[n0]-=3;
    if (k=='n' || k=='m') {
    }

    if (k=='0') reset_vars();
}

void reshape(int x, int y) {
    width = x;
    height = y;
    initPixelBuffer();
    glViewport(0, 0, x, y);
    glLoadIdentity();
    glOrtho(0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
}

void cleanup() {
    cudaGLUnregisterBufferObject(pbo);
    glDeleteBuffersARB(1, &pbo);
}


int main( int argc, char** argv) {
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
    glutInitWindowSize(width, height);
    glutCreateWindow("CUDA WhiteNoise");
    glutDisplayFunc(display);
    glutKeyboardFunc(keyboard);
    glutReshapeFunc(reshape);

    glewInit();
    if (!glewIsSupported("GL_VERSION_2_0 GL_ARB_pixel_buffer_object")) {
        fprintf(stderr, "OpenGL requirements not fulfilled !!!\n");
        exit(-1);
    }
    
    reset_vars();
    
    makeRasterFont();
    initPixelBuffer();

    atexit(cleanup);
    glutMainLoop();
    return 0;
}

void reset_vars() {
    h_vars[a0] = -0.4;
    h_vars[a1] = -1.0;
    h_vars[a2] = -0.4;
    h_vars[a3] =  0.1;

    h_vars[b0] =  0.76;
    h_vars[b1] = -0.4;
    h_vars[b2] =  0.4;
    h_vars[b3] =  0.76;

    h_vars[fw] = 1.0f;
    h_vars[fh] = 1.0f;

    h_vars[n0] = 80;

    keyb_step = 0.02;
}
